A0-A1S8  977 


UNCLASSIFIED 


A  FAST  ALGORITHM  FOR  NON-NEHTONIAN  FLOM(U>  WISCONSIN 
UNIV-HADISON  MATHEMATICS  RESEARCH  CENTER  D  S  NALKUS 
SEP  85  HRC-TSR-2868  DAAG29-88-C-A841 


m 


UNIVERSITY  OF  WISCONSIN 
MATHEMATICS  RESEARCH  CENTER 

A  FAST  ALGORITHM  FOR  NON-NEWTONIAN  FLOW 


David  S.  Malkus 

Technical  Summary  Report  ^2860 
September  1985 

ABSTRACT 


The  goals  of  the  project  described  here  are  twofold:  First,  to  turn  an  existing  pilot  al¬ 
gorithm  for  the  steady  flow  of  non-Newtonian  memory  fluids  into  a  robust  and  efficient 
algorithm.  Second,  render  enhancements  of  the  method’s  current  capabilities  computa¬ 
tionally  feasible.  Such  enhancements  include  fully  coupled  thermal  dependence,  material 
compressiblity,  and  free  surface  flows.  The  pilot  algorithm  is  a  finite  element  method  whose 
novelty  lies  in  its  computation  of  the  stress  field  in  a  nonlinear  iteration  scheme.  The  stress 
at  a  point  is  a  non-local  functional  of  the  current  velocity  iterate,  and  the  pilot  method 
has  demonstrated  the  feasiblity  of  reliable  computation  with  such  constitutive  equations. 
Before  the  method  can  take  its  place  as  a  reliable  scientific  and  engineering  tool,  intensive 
effort  must  be  made  to  reduce  the  computational  cost  in  the  manner  described  here. 
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A  FAST  ALGORITHM  FOR  NON-NEWTONIAN  FLOW 

David  S.  Malkus 


1.  VISCOELASTIC  FLUIDS.  The  following  equations  are  solved  numerically,  using  the 
finite  element  method  .1  —  4  :  The  equat  ions  of  steady  motion. 


V  a 


p( u  •  V)u 


where  u  is  the  velocity  field,  a  the  stress  tensor,  f  a  body  force,  and  p  the  density.  The 
equation  of  continuity  for  an  incompressible  fluid  is 

^  •  u  -  0  (2) 

For  non-Newtonian  fluids,  the  crucial  equation  is  the  constitutive  equation, 

a  -  -pi  +  2p(0)/?e  +  (1  R)o'  (3) 

where  p  is  an  isotropic  contribution  to  the  stress,  p(0)  is  a  zero-shear  viscosity,  R  is  a  ratio 
of  a  retardation  time.  A,  to  a  retardation  time,  T ,  and  a '  is  an  extra  stress  tensor.  The 
ratio,  R,  and  its  complement  determine  the  proportion  of  the  stress  which  is  Newtonian 
—  and  usually  is  the  result  of  a  Newtonian  solvent  —  and  the  complementary  proportion 
from  the  extra  stress  —  usually  due  to  long-chain  molecules  (such  as  polymers)  disolved 
in  the  Newtonian  solvent. 

There  are  many  proposed  forms  for  the  extra-stress  tensor;  there  are  two  basic  cat¬ 
egories:  the  differential  and  integral  models  'ij.  Here  we  shall  only  be  concerned  with  the 
integral  form. 
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where  pf,  is  a  constant  determining  p(0),  and  is  a  strain  measure,  measuring  the 

deformation  which  carried  the  particle  from  its  position  at  time  r  in  the  past  to  the  stress 
evaluation  point  at  the  present  time.  0.  The  strain  measures  are  the  same  kind  employed 
in  finite  elasticity.  The  memory  functions,  m/,  are  usually  sums  of  exponentials,  each  with 
amplitude  determined  by  the  modulus,  C and  decay  constant,  pk ,  which  determines  the 
fraction  of  the  basic  decay  rate,  T. 
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Thus  a  computational  method  must  determine  the  deformation  history  of  every  stress  eval¬ 
uation  point  required  to  solve  the  equations  of  motion  in  some  approximate  way,  compute 
the  required  strain  measure  —  which  is  almost  always  highly  nonlinear  in  its  dependence 
on  the  velocity  field  —  and  then  approximate  the  history  integral  over  an  infinite  interval. 
This  just  computes  the  stress,  and  then  the  stress  computation  must  be  imbedded  in  some 
iterative  scheme  to  produce  an  approximate  solution  to  the  highly  nonlinear  equations  of 
mot  ion . 

11.  SOME  JFLOWS.  In  this  section  we  give  a  brief  description  of  some  of  the  flows  to  which 
the  current  method  is  being  applied.  The  geometry  of  these  flows  is  quite  simple  and  the 
results  obtained  do  not  illustrate  the  real  power  of  the  finite  element  method.  It  is  hoped 
that  the  reader  will  appreciate  that  the  method  described  here  is  still  very  much  in  the 
development  stage,  and  that  the  problems  so  far  investigated  by  the  author  and  other 
researchers  are  intended  to  isolate  the  complexity  inherent  in  the  non-Newtonian  nature 
of  the  flow  from  other  possible  complications.  Nevertheless,  there  seems  to  be  a  good  deal 
of  physical  interest  in  the  problems  pictured  here,  in  spite  of  their  geometric  simplicity. 


FIGl'RE  1 

f  low  over  a  transverse  slot  computed  by  the  pilot  method  on  a  mesh  of  1008  elements. 

1  he  first  flow  is  a  plane  flow  over  a  transverse  slot.  The  streamlines  plotted  in  Figure 
1  are  taken  from  a  solution  computed  b>  the  author,  using  a  constitutive  equation  of  his 
own  devising  |  and  the  mesh  of  1008  crosed-triangle  macroelements  illustrated  in  ref.  1. 
I  he  flow  is  at  a  "Deborah  number  of  1.7  (this  can  be  thought  of  as  a  dimensionless  shear 
rate),  flow  is  from  right  to  left,  and  undisturbed  flow  profiles  have  been  imposed  at  the 


inflow  and  outflow.  Actually,  since  there  is  fluid  memory,  the  inlet  condition  is  that  the 
flow  continues  forever  upstream  as  undisturbed  plane  Poiseuille  flow.  Figure  1  illustrates  a 
chracteristic  tilt  to  the  vortex  in  the  slot,  which  is  opposite  in  direction  to  the  tilt  observed 
in  Newtonian  flows  with  non-zero  Reynolds  number  3  . 

The  interest  in  flows  over  transverse  slots  aris'  from  the  fact  that  there  seems  to  be 
an  important  relation  between  the  difference  between  the  pressures  at  top  and  bottom  of 
the  slot  and  the  first  normal-stress  difference  of  the  fluid  in  the  undisturbed  flow  1  -  3  . 
1  here  seems  to  be  a  discrepancy  between  what  the  numerical  models  predict  and  laboratory 
experiments  measure  in  such  flows,  and  it  is  one  of  the  author  s  highest  priorities  to  resolve 
this.  The  results  could  have  important  ramifications  for  devices  designed  to  measure  the 
first  normal-stress  difference  using  “hole-pressure”  measurements. 


FIGURE  2 

Abrupt  contraction  flow  computed  by  the  pilot  method  with  940  elements. 

The  next  flow  is  that  of  flow  through  an  abrupt,  planar  contraction.  Figure  2  pictures 
such  a  flow:  flow  is  from  left  to  right  so  that  the  fluid  is  being  forced  form  the  larger  to 
the  smaller  channel.  Because  symmetry  is  assumed,  the  computational  domian  is  only  the 
top  half  of  a  channel  cross-section.  The  flow  pictured  here  uses  the  same  fluid  model  as 
that  of  figure  1.  at  a  slightly  lower  Deborah  number.  3.7.  Inflow  and  outflow  conditions 
are  imposed  as  before;  extreme  care  is  taken  to  match  the  flux  at  inflow  and  outflow 
boundaries. 

The  interest  in  this  flow  stems  from  the  fact  that  some  fluids  seem  to  behave  quite 
different!)  than  others  in  contraction  flow.  Some  fluids,  such  as  polystyrene  or  high-density 


polyethelene  melts,  seem  to  have  relatively  smaller  “dead-spaces"  or  recirculation  regions 
at  high  Deborah  number  than  at  low  Deborah  number,  while  some  branched  polymers, 
such  as  low-density  polvethelyne.  seem  to  do  quite  the  opposite,  developing  recirculation 
regions  emanating  from  entry  which  dominate  the  whole  flow-field.  The  flow  pictured  in 
Figure  2  has  about  the  same  size  recirculation  region  as  a  flow  of  the  same  fluid  at  low 
shear  rate.  The  author  is  interested  in  further  study  of  this  flow  in  order  to  find  out  what 
property  of  the  constitutive  mode!  is  associated  with  entry  vortex  behavior.  The  ability 
to  predict  recirculation  size  is  of  practical  import  because  fluid  trapped  inside  dead-spaces 
tends  to  degrade.  It  would  be  useful  to  able  to  determine  how  much  polymer  degradation 
could  be  expected  in  a  given  die  as  a  function  of  measurable  material  properties. 
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FIGURE  3 

Newtonian  flow  in  a  plane  cross-section  of  a  journal-bearing. 

The  final  flow  is  one  for  which  the  author  has  as  yet  no  results;  that  is  the  flow  in 
a  journal-bearing.  Figure  pictures  the  cross-section  of  two  eccentrically  placed  cylinders 
with  fluid  between  them  to  lubricate  and  prevent  solid-to-solid  contact.  There  are  two 
important  aspects  to  this  flow  which  differ  from  the  previous  two  flows;  First,  there  are 
no  inflows  and  outflows,  and  second,  there  are  no  domain  corners  to  generate  singularities 
in  t lie  st ress- field. 

The  author’s  interest  in  this  problem  is.  first,  just  that  it  is  quite  different  from  the 
other  two.  It  will  be  interesting  to  observe  the  behavior  of  the  numerical  method  here 
because  it  omits  two  puzzling  aspects  of  memory-fluid  problems:  history-dependence  at 
inlets,  and  stress-singularities  of  unknown  character.  There  is  also  physical  interest  in  this 
problem  because  the  effects  of  fluid  elasticity  on  the  load-bearing  capacity  of  the  bearing 


may  be  beneficial.  It  would  be  able  to  predict  load  bearing-capacity  from  measurable 
materia!  properties,  and  numerical  modelling  may  help  to  do  so. 

Ill.  COMPUTATIONAL  METHOD.  The  strain  measure  in  the  integrand  of  eq.  (1)  is 
assumed  to  be  determined  by  a  deformation  gradient.  Eo(r).  just  as  in  nonlinear  elasticity. 
Only  in  the  present  case,  the  deformation  gradient  is  assumed  to  be  computable  from  a 
system  of  linear,  non-constant  coefficient  ordinary  differential  equations  along  the  path 
followed  by  each  particle  at  which  the  stress  is  to  be  evaluated.  The  usual  deformation 
gradients  of  large-strain  elasticity  can  also  be  obtained  from  special  cases  of  the  following 
evolution  equations: 

x(r)  =  v  \{r) 
x(0)  -  X0 

(5) 

E0(')  —  F(x(r).Vv  x(r);)E0(r) 

Eo(0)  -  I 

The  first  two  sets  of  equations  determine  the  pathline  (streamline)  followed  by  a  particle 
to  bring  it  from  its  position,  x.  at  time  r  in  the  past  to  its  present  position  at  the  stress 
evaluation  point,  Xo-  To  evaluate  the  integrand  of  eq.  (4),  these  equations  are  solved  as 
an  initial  value  problem  in  reverse  time.  This  determines  the  non-constant  coefficient  in 
the  evolution  equation  for  the  gradient,  which  is  assumed  to  be  a  traceless  matrix,  F.  The 
common  deformation  gradient.  is  obtained  when  F  is  Vv  itself. 

The  fundamental  strategy  of  the  current  numerical  method  is  to  choose  constant 
strain-rate  finite  elements:  then  the  evolution  equation  is  a  const  ant -coefficient  equation 
on  each  element.  This  strategy  is  enabled  by  a  basic  property  of  linear  ODEs:  If  we  define 
a  deformation  gradient.  Er,.  relative  to  time  r j  by 

Er,  =  FEt, 

Er,(r,)  =  I 

then  the  strain  relative  to  the  present  time,  evaluated  at  any  earlier  time  is  give  by  matrix 
multiplication: 

E0(r)  =  E0(r,)ETl(r)  (7) 

This  provides  interface  conditions  between  finite  elements,  so  that  only  constant-coefficient 
equations  need  be  solved  on  each  element:  it  turns  out  that  such  solutions  are  known 
analytically,  as  is  the  pathline  and  transit  time  along  it  1  —  3  . 

Thus,  given  an  estimate  of  the  solution  to  the  problem  it  terms  of  a  velocity  field, 
the  integrand  of  eq.  (4)  can  straight-forwardly  be  computed  at  each  historical  time.  In  the 
current  method,  this  is  used  in  conjunction  with  a  specially  devised  Gaussian  quadrature 
formulas  to  approximate  o': 


(6) 


With  what  we  have  thus  far.  the  stress  can  be  computed  in  any  trial  velocity  field;  to 
approximate  the  solution  to  eqs.  (1)  and  (2).  the  usual  Galerkin  procedure  can  be  followed 
in  which  the  residual  of  eq.  (1)  is  doited  into  a  test  function.  \h ,  drawn  from  the  same 
space  as  the  trial  solutions,  and  the  result  is  integrated  over  the  problem  domain.  After 
integration  by  parts  and  replacement  of  the  spatial  integral  by  a  numerical  integral  with 
points  and  weights  6f.  we  get  something  which  looks  like 

]T0t  o'  ■  Vv*  -  2z(V  i/)(V  •  v'  )  -  p( u  •  V)u''  ■  vh  -  vh  f  (&)  =  0  (9) 

The  pressure  term  ofeq.  (3)  has  been  replaced  by  a  penalty  term  [3  with  penalty  parameter 
z\  thus  there  are  no  explicit  pressure  unknowns,  and  ti.e  continuity  equation  (2)  is  satisfied 
to  0{z  ').  Eq.  (9)  illustrates  the  R  -  0  case;  for  nonzero  R.  the  obvious  modification  of 
adding  a  Newtonian  viscous  term  is  made. 

The  important  point  to  observe  about  eq.  (9)  is  that  to  evaluate  its  residual,  it  is 
required  to  evaluate  the  stress  at  the  points  (,  by  means  already  discussed.  To  complete  the 
method,  what  is  needed  is  a  means  of  correcting  estimates  of  the  discrete  solution,  based  on 
evaluation  of  the  residual:  Newton's  method  might  an  example  of  such  a  procedure,  but,  as 
we  shall  see.  this  is  not  entirely  straight-forward.  The  current  algorithm  employs  the  inverse 
Brovden  method  1.2  to  solve  the  discrete  nonlinear  equations.  An  important  point  to  be 
made  her  '  is  that,  regardless  of  the  choice  of  iterative  scheme,  the  method  outlined  here  is 
enormously  costly  in  practice,  because  for  a  reasonably  fine  mesh,  each  evaluation  of  the 
stress-field  values  at  the  spatial  integration  points  is  a  potentially  formidable  computation. 

IV.  FAST  ALGORITHMS..  The  method  outlined  in  the  previous  section  applies  to  isother¬ 
mal.  incompressible  flows  in  a  fixed  spatial  domain.  These  restrictions  are  not  essential; 
material  compressibility  and  temperature  dependence  can  be  handled  in  very  similar  fash¬ 
ion  if  "artificial  (historical)  time1'  is  introduced,  in  which  either  density  or  temperature 
are  used  to  change  the  time  variable  along  the  pathlines  in  such  a  way  that  a  traceless 
matrix  in  the  evolution  equation  is  obtained  5  .  The  transformation  to  artificial  time  does 
not  in  itself  seem  to  be  computationally  costly,  but  these  problems  involve  added  levels  of 
complexity  to  an  already  complicated  solution  procedure  with  additional  fields  and  corre¬ 
sponding  equations.  The  resulting  phenomena  are  likely  to  be  more  intricate  in  detail  and 
more  nonlinear  in  character.  A  similar  observation  can  be  made  about  free-surface  flows: 
a  well-developed  methodology  exists  6  to  solve  such  problems,  which  can  be  directly  in¬ 
terfaced  with  the  method  outlined  here,  but  this  also  certain  to  render  the  computations 
more  formidable  than  they  now  are.  With  the  current  algorithm,  computations  on  a  mesh 
which  is  refined  only  to  the  extent  which  seems  to  be  required  to  obtain  acceptable  accu¬ 
racy.  at  a  shear  rate  normally  occuring  in  polymer  processing,  the  computation  of  a  steady 
solution  can  take  as  long  a  40  minutes  on  a  (’ray  1-A.  This  must  be  reduced  drastically  if 
the  algorithm  is  to  used  routinely  in  scientific  and  engineering  research,  particularly  if  the 
more  physically  realistic  enhancements  mentioned  above  are  to  be  added.  The  remainder 
of  this  paper  will  discuss  several  approaches  to  the  reduction  in  computational  cost  which 
are  currently  being  implemented  or  investigated  by  the  author. 


Vectorization  of  Linear  Equation  Solving.  A  variety  of  new  computers  have  the  capability 
of  carrying  out  hardware  vector  operations:  rearranging  the  computer  code  in  such  a  way 
that  the  compiler  can  take  advantage  of  this  capability  can  result  in  substantial  savings  in 
computational  cost.  One  part  of  a  typical  code  where  such  savings  have  a  good  chance  of 
being  realized  is  in  the  solution  of  linear  equations.  Unfortunately,  in  the  current  algorithm, 
it  is  not  expected  that  this  can  dramatically  reduce  the  run  time.  Linear  equations  are 
solved  in  the  nonlinear  iteration  scheme,  but  this  appears  to  account  for  a  small  portion  of 
the  computational  cost.  The  major  portion  of  the  calculation  is  carried  out  at  the  element 
level  with  small  arrays  or  scalar  quantities  involved  in  resolving  element  boundary  crossings 
and  accumulating  the  deformation  gradient  by  small  matrix  multiplication.  Vectorization 
offers  little  hope  of  speeding  up  these  calculations.  On  the  the  other  hand,  it  is  expected 
that  linear  equation  solving  will  begin  to  play  a  more  and  more  important  role  with  the 
planned  enhancements  to  the  code  discussed  earlier.  The  Jacobian  terms  corresponding  to 
the  thermal  energy  equation  are  easy  to  form;  likewise  the  part  of  the  Jacobian  associated 
with  inertial  terms  and  the  unknowm  free-surface  transformation  are  easy  to  deduce.  Also, 
active  research  is  under  way  aimed  at  producing  the  Jacobian  terms  associated  with  the 
non-Newtonian  viscous  terms  (see  below). 

In  short,  the  future  development  of  the  code  seems  to  point  in  the  direction  of  an 
algorithm  which  has  a  large,  unsymmetric.  and  possibly  not  banded  matrix  to  factor  at 
each  one  of  dozens  of  possible  iterates.  The  current  iteration  method  has  only  one,  banded, 
symmetric,  positive-definite  matrix  to  factor  at  the  outset,  and  a  back-substitution  at  each 
iteration  (the  unsymmetric  Jacobian  contribution  of  the  inertial  terms  is  left  to  the  inverse 
updating  scheme).  It  therefore  seems  appropriate  to  modify  the  code  at  the  present  time 
to  take  full  advantage  of  vectorization.  in  order  to  make  sure  the  linear  equation  solving 
phase  remains  in  the  bacground,  as  it  should. 

Adaptive  Memory  Quadrature.  The  area  which  seems  to  show  most  promise  in  reduction  of 
the  the  computational  cost  is  that  of  the  stress  calculation  at  an  individual  stress  evaluation 
point .  There  seem  to  be  several  possible  approaches,  the  underlying  strategy  of  all  of  them 
is  to  take  advantage  of  the  fact  that  the  stresses  are  being  evaluated  in  what  is  hoped  will 
be  a  comergent  sequence  of  velocity  iterates.  Particuarly  further  along  in  the  sequence, 
the  previous  iterate  should  be  able  to  provide  a  guide  to  estimate  how  much  computation 
is  absolutely  necessary  at  the  next  iteration. 

Perhaps  the  most  obvious  way  to  do  this  is  to  use  the  previous  iterate  to  determine 
what  ,\;  of  eq.  (H)  should  be  in  the  next  iterate.  It  is  observed  that  in  some  flows,  very 
many  fewer  quadrature  points  are  needed  to  accurately  compute  o’  than  in  other  flows. 
1  he  strategy  will  be  to  begin  the  iterations  with  a  nominal  number  of  points  for  each  stress 
evaluation  point  and  increase  or  decrease  that  number  in  succeeding  iterations,  based  on 
adaptive  criterea  determined  from  previous  iterations. 

Jacobian  Approximation.  The  reason  that  the  present  algorithm  employs  an  updating 
scheme  rather  than  a  direct  calculation  of  the  Jacobian  is  that  it  is  not  a  trivial  mat¬ 
ter  to  construct  the  Jacobian,  or  even  write  down  a  closed  form  expression  for  it.  It  is 
clear  that  the  stress  at  a  point  can  depend  on  velocities  far  from  that  point,  and  therefore 


the  Jacobian  of  the  residual  of  eq.  (9)  cannot  have  the  usual  finite  element  band  and  or 
sparsity  structure.  In  ref.  7.  an  approximation  scheme  for  the  Jacobian  is  proposed.  It  is 
not  clear  at  present  whether  this  approximation,  some  other,  or  even  an  exact  computation 
of  the  Jacobian  is  best  (the  latter  may  be  possible  to  undertake  —  it  is  not  clear  at  this 
time).  But  the  work  of  ref.  7  shows  clearly  the  complexity  involved.  The  terms  of  the 
Jacobian  contribution  from  the  extra  stress  are  computed  by  tracking  along  streamlines. 
The  resulting  Jacobian  element  matrices  are  noi  square:  Their  column  dimension  depends 
on  the  number  of  different  elements  the  particle  path  passes  through  before  the  final  in¬ 
tegration  point  of  eq.  (8)  is  located.  It  appears  that  a  frontal  solution  technique  is  called 
for  in  order  lo  handle  the  resulting  global  matrix  7  . 

INeudo-Dynarnic  Relaxation.  The  hope  in  computing  the  Jacobian  is  that  the  computa¬ 
tional  cost  will  be  more  than  returned  in  improved  convergence  rate  over  the  inverse  Brov- 
den  algorithm  afforded  b\  Newton's  method  or  modified  Newton's  method  with  Broyden 
updates.  But  t  fie  complexity  of  the  Jacobian  calculation  is  such  that  this  may  never  be  real¬ 
ized.  and  it  is  well  worth  the  investigation  of  other  improvements  to  the  iterative  solution  of 
the  nonlinear  equations.  One  avenue  currently  being  explored  is  that  of  “pseudo-dynamic 
relaxt ion."  The  problem  is  cast  as  a  time-dependent  problem  and  steady  solutions  are 
obtained  b>  letting  the  transient  phenomena  die  out.  In  the  algorithm  presented  here, 
the  transient  behavior  does  not  represent  the  true  dynamics  of  the  non-Newtonian  fluid; 
the  stress  is  computed  in  the  current  velocity  field  as  if  it  had  been  a  steady  field  for 
all  time,  hence  the  name  "pseudo-dynamic."  To  do  otherwise  would  involve  complexities 
beyond  what  seems  manageable  at  present,  though  implementation  of  the  pseudo-dynamic 
algorithm  does  open  the  door  for  future  exploration  of  true  dynamic  behavior. 

One  may  easily  verify  that  the  steady-states  of  the  pseudo-dynamic  algorithm  are 
the  same  as  the  steady  states  of  the  true  dynamic  algorithm.  The  reason  for  taking  the 
pseudo-dynamic  approach  is  to  produce  a  different  kind  of  steady-state  iteration  scheme, 
in  which  the  damping  of  the  high  frequency  modes  in  the  pseudo-dynamic  response  can 
be  controlled  by  choice  of  time-stepping  method  and  pseudo-lime  step.  The  reason  that 
this  seems  to  be  a  worthwhile  avenue  to  explore  is  suggested  by  recent  work  of  \ .  Renardv 
and  M.  Renardv  8  .  They  found  that  with  a  certain  spatial  discretization  of  the  linearized 
operator  associated  with  the  equations  of  motion  of  a  Maxwell  fluid  in  a  shearing  flow, 
there  were  apparently  spurious  eigenvalues  extremely  close  to  the  right  half-plane,  evi¬ 
dently  introduced  by  the  discretization.  If  this  were  also  a  consequence  of  finite  element 
discretization,  there  could  be  severe  consequences  for  iterative  methods  which  behave  like 
temporal  iteration  schemes.  It  is  hoped  that  by  controlling  the  time  step  and  parameters 
of  the  pseudo-dynamic  time-stepping  method,  the  damping  of  the  high  frequency  modes 
associated  with  any  spurious  eigenvalues  can  be  damped  to  produce  nearer  monotonic, 
more  rapid  convergence  of  the  resulting  iterative  method.  There  is  no  worry  hereof  damp¬ 
ing  out  interesting  transient  behavior  the  transient  behavior  is  not  correct,  and  all  that 
it  is  required  is  that  it  be  damped  out  as  rapidly  as  possible. 

The  following  algorithm  is  based  on  Hughes.  Liu.  and  Brooks  pedictor-corrector 
algorithm  for  the  Navier-Stokes  equat ions  9  .  but  with  the  possibility  of  more  fully  implicit 


inner  iterations  at  each  time  step: 


Man.i  -  Cvn*i  -  N(vn^j)  -  Q(vn.  x)  -  Fn.x 

Q  /  B TodV  -  Cv 

Jci 

N  nonlinear  inertia!  term.  exc).  time  deriv. 
vn°)i  vn  -  (1  -  ~i)&tan 

vn*iJ  "  Vn-1  -  +  TfA/CJv^!  - 

+  '>AfN(v^1)-Mv^)1  +  ')A<Fn.1} 

an-  1  (V„.  x  V^Jj)/*) At 

J  M  -  (-vAf^)  -r 

'  A\r  '  Opt  '  Q'y  '  opt 


(10) 


M  is  the  finite  element  “mass  matrix.”  C  the  Newtonian  viscous  and  penalty-pressure 
matrix,  and  Fn_i  the  applied  force  vector  at  time  step  n  +  1.  B  is  the  usual  finite 
element  matrix  of  shape  function  derivatives  and  a  is  the  stress,  computed  in  as 

described  in  previous  sections:  vn^.i  without  the  superscript  of  inner  iteration  is  the  “fully 
converged”  result  of  inner  iteration  at  time-level  n-1.  Choice  of  the  number  of  inner 
iterations  is  open,  so  that  vn_x  could  result  from  just  one  correction  cycle,  or  many.  An 
important  aspect  of  eq.  (10)  is  found  in  those  terms  labelled  by  (-)opt;  a  fully  implicit 
treatment  would  employ  exact  Jacobian  terms  here.  At  the  other  extreme  is  Hughes,  Liu, 
and  Brooks  method:  they  use  C  to  approximate  both  of  these  terms  and  do  only  one  inner 
iteration.  The  present  non-Newtonian  implementation  uses  C  initially,  updated  by  the 
inverse  Broyden  method  during  a  number  of  inner  iterations.  If  it  proves  to  be  effective, 
Newton  or  modified  Newton  Broyden  iterations  could  be  employed  in  the  inner  iterations. 
It  is  instructive  to  note  that  the  direct  steady-state  Broyden  algorithm  mentioned  earlier 
is  obtained  as  a  special  case  of  eq.  (10)  with  -y  -  1  and  an  infinite  time  step. 

V.  SUMMARY.  A  pilot  numerical  method  for  the  computation  of  solution  to  memory 
fluid  flow  problems  has  been  described.  This  method  has  shown  that  such  computations 
are  feasible  but  extremely  costly.  More  reasonable  physical  assumptions  than  those  of 
isothermal,  incompressible  flow  in  a  fixed  domain  are  on  the  near  horizon  but  are  bound 
to  increase  the  computational  cost.  A  number  of  ways  of  improving  the  computational 
performance  of  the  algorithm  have  been  proposed  here  and  are  in  the  implementation 
stage.  These  improvements  will  go  together  to  make  what  the  author  refers  to  as  “a  fast 
algorithm  for  non-Newtonian  flow.”  It  is  hoped  that  this  fast  algorithm  can  transform  the 
method  described  here  from  pilot  code  to  useful  computational  tool  for  ihe  investigation 
of  problems  in  viscoelasticity  and  rheology. 
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